home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / Quick3 / Quick3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-04  |  14.0 KB  |  349 lines  |  [TEXT/KAHL]
  1. /*
  2. Quick3.c 
  3. Copyright (c) 1990-1993 Denis G. Pelli 
  4. A C program to make maximum likelihood fits of simple models to
  5. psychometric data. This is a replacement for, and to some extent is
  6. derived from the old QUICK.FOR program written by A.B.Watson. See
  7. A.B.Watson (1979) Probability summation over time. Vision Research 19, 515-522.
  8.  
  9. The principal difference between this program and the original QUICK, is in how
  10. the goodness of fit is calculated. Both programs assess the goodness of fit
  11. by comparing the log likelihoods of the model fit (usually Weibull) with
  12. a null hypothesis. Minus two times this log likelihood ratio is approximately
  13. chi square, with a number of degrees of freedom equal to the difference in 
  14. degrees of freedom of the two hypotheses. QUICK used an unconstrained psychometric
  15. function, which adopted the actual proportion correct at each contrast. Quick3 uses
  16. a monotone hypothesis which is constrained to be monotonically increasing as a function
  17. of contrast. The virtue of the latter hypothesis is that it is easier to estimate
  18. its degrees of freedom when there are unequal numbers of data points at each contrast.
  19. (However, this estimate is seat-of-my-pants.) Finally, this program actually computes
  20. the significance for you, from the chi square value and the estimate of its degrees of
  21. freedom. This significance value is of some use in choosing a small fraction of your
  22. fits to discard (e.g. if significance is less than 0.05). This significance is
  23. misleadingly low when both hypotheses have nearly the same number of degrees of
  24. freedom (i.e. you only have a few contrasts). Looking at the significance is 
  25. helpful, but not a substitute for thinking about whether your results are really
  26. reasonable. 
  27.  
  28. Note that Quick3.c is just a front end, dealing with the user interface and
  29. reading and writing files. All the work is done by one call to PsychometricFit().
  30. In many cases it will be convenient to call that routine directly from 
  31. your program that collects (or generates) the data, rather than saving the data
  32. in a file for subsequent analysis by Quick3.c. If you want to use a psychometric
  33. model other than the Weibull function all you need to do is write a small function
  34. to implement it, based on Weibull.c. It should be easy. PsychometricFit()
  35. uses whatever psychometric function is supplied in the call.
  36.  
  37. Quick3 reads in your data from a text file and analyzes them. 
  38. The results are shown on the screen, and are saved in two kinds of
  39. output file. The .fit file is in Excel format and has just a minimal one-line
  40. summary for each condition, giving the parameters and goodness of fit. The .plot
  41. file is in CricketGraph format and is suitable for plotting the psychometric
  42. data along with the Weibull and Monotonic fits. For Macintosh users, the output files
  43. are TEXT files, but have the appropriate creator so that double-clicking them
  44. opens the data file in the appropriate application, either Excel or Cricket Graph.
  45.  
  46. You should also have received a CricketGraph format file called "CricketGraph Quick3".
  47. Put this file in the same folder as the CricketGraph application. When CricketGraph
  48. starts up it will read in the format file, and later you'll be able to select that
  49. format from CricketGraph's Format menu when you plot your data.
  50.  
  51. Quick3's input file format is simple and flexible. It is a text file. 
  52. Lines beginning with "#" are treated as comments. Comments are either copied directly to
  53. the .fit file, or used as a name for the next "condition" (i.e. block
  54. of data). The file must begin with at least one comment line. The last of several
  55. contiguous comment lines will be taken to be the name for the following condition.
  56. All the non-comment lines until the next comment will be interpreted as the data
  57. for that condition.
  58.  
  59. An experiment usually consists of several conditions, or blocks of trials.
  60. Within the condition the observer will be tested repeatedly at various contrasts.
  61. The results consist of the number of trials correct at each contrast. The results
  62. may be described by a set of "contrast records". Each contrast record
  63. consists of a contrast (e.g. 0.012), the number of trials at that contrast (e.g. 1 or
  64. 100), and the number of correct responses (e.g. 0 or 79). Quick3 expects a contrast
  65. record to look like this:
  66. 0.012 100 79
  67. The three numbers are all on one line, and separated by white space,
  68. so sscanf can parse them. Anything on the line after the third number is ignored.
  69. All contrast records until the next comment (line beginning with "#") or the end of file
  70. are assumed to be from the same condition and are analyzed together. The contrast 
  71. records may be in any order. The records will be sorted into order of ascending contrast.
  72. You may have any number of trials at each contrast, up to about two billion (i.e. 
  73. LONG_MAX in limits.h). It's okay to have multiple records with equal contrasts. The data
  74. will be merged, adding up the trials at that contrast. You can even be really lazy and
  75. just write out each trial as a separate contrast record, in whatever order, without
  76. keeping track of how many, etc. 
  77.  
  78. You may add records with zero trials (and zero correct) so as to have fitted
  79. values computed at those contrasts in the .plot file.
  80.  
  81. Here's an example of a data file, suitable for analysis by Quick3:
  82. #This is a sample file. 
  83. #The first condition is called "monocular" and has 60 trials at 6 contrasts.
  84. #monocular
  85. 0.1 10 4
  86. 0.2 10 5
  87. 0.3 10 7
  88. 0.4 10 6 Anything after the specified data is ignored.
  89. 0.5 10 8
  90. 0.6 10 10
  91. # The "binocular" condition is bigger: 1000 trials.
  92. #binocular
  93.  0.056     100    53
  94.  0.064     100    53
  95.  0.073     100    65
  96.  0.083     100    77
  97.  0.094     100    81
  98.  0.107     100    84
  99.  0.121     100    96
  100.  0.138     100    99
  101.  0.156     100   100
  102.  0.178     100   100
  103. #all done
  104.  
  105. Here's the resulting .fit file (tabs will space properly in Excel, but not here):
  106. Condition    logAlpha    beta    gamma    delta    free params    signif.    Chi sq.    d.f.    trials    contrasts
  107. #This is a sample file. 
  108. #The first condition is called "monocular" and contains 60 trials at 6 contrasts.
  109. monocular    -0.305    8.269    0.524    0    4    0.183    1.773    1    60    6
  110. # The "binocular" condition is bigger: 1000 trials.
  111. binocular    -1.015    3.811    0.482    0    4    0.479    5.521    6    1000    10
  112. #all done                                        
  113.  
  114. SOURCES:
  115.  
  116. Quick3.h
  117. LogLikelihood.c
  118. MonotonicFit.c
  119. PsychometricFit.c
  120. Quick3.c
  121. SortAndMergeContrasts.c
  122. Weibull.c
  123. #From Denis Pelli's VideoToolbox:
  124. VideoToolbox.h
  125. Binomial.c
  126. ChiSquare.c
  127. Normal.c
  128. SetFileInfo.c        # Used only on the Macintosh
  129. #From Numerical Recipes in C:
  130. nr.h
  131. NRUTIL.h
  132. BRENT.C
  133. F1DIM.C
  134. LINMIN.C
  135. MNBRAK.C
  136. NRUTIL.C
  137. POWELL.C
  138.  
  139. HISTORY:
  140. 4/7/90        dgp    wrote it.
  141. 4/10/90        dgp    changed .fit format from %.3f to %.4f. Accept only printing characters
  142.                 in the condition name, ignoring everything after the first non-printing
  143.                 character. This deals with the fact that if an Excel file is used for
  144.                 input, there are lots of trailing tabs that should be ignored. Checked
  145.                 that files were actually open before closing 'em. Jeesh!
  146. 10/29/90    dgp    tidied up comments.
  147. 1/20/91        dgp    updated calls to BinomialUpperBound & BinomialUpperBound to
  148.                 conform to new definition.
  149. 8/24/91        dgp    Made compatible with THINK C 5.0.
  150. 11/17/92    dgp Now SetVol() after calling SFGetFile, so that file can be in any folder,
  151.                 not just the same folder that Quick3 is in.
  152. 1/19/93        dgp    put #ifs around the Mac-dependent code.
  153. 2/20/93        dgp    added call to Require().
  154. */
  155. #include "VideoToolbox.h"
  156. #include "Quick3.h"
  157. #include <nr.h>                /* Numerical Recipes in C*/
  158. #if MACINTOSH
  159.     #include <QuickDraw.h>
  160.     #include <StandardFile.h>
  161. #endif
  162.  
  163. void main(void)
  164. {
  165.     static dataRecord data,monotonicData;
  166.     contrastRecord *cPtr;
  167.     static paramRecord params, initParams;
  168.     double *paramPtr=¶ms.logAlpha;    /* a generic way of accessing the parameters */
  169.     int i;
  170.     double chiSquare,modelLL,monotonicLL;    /* log likelihood */
  171.     int chiSquareDF,modelDF,monotonicDF;    /* degrees of freedom */
  172.     double significance;
  173.     FILE *dataFile=NULL,*fitFile=NULL,*plotFile=NULL;
  174.     static char dataFileName[50],fitFileName[50],plotFileName[50],string[100];
  175.     static char conditionName[50];
  176.     long trials;
  177.     unsigned int a;
  178.     PsychometricFunctionPtr ModelFunction=&Weibull;    /* a psychometric function */
  179.     static const char modelName[]="Weibull";
  180.     static const char paramName[PARAMS][16]={"logAlpha","beta","gamma","delta"};
  181.     #if MACINTOSH
  182.         Point where;
  183.         static SFTypeList typeList;
  184.         static SFReply reply;
  185.     #endif
  186.  
  187.     /* initial values for the parameters of the psychometric function */
  188.     params.logAlpha=0.0;
  189.     params.beta=3.5;
  190.     params.gamma=0.5;
  191.     params.delta=0.01;
  192.     printf("\n");    /* ask THINK C to init QuickDraw */
  193.     
  194.     #if MACINTOSH
  195.         Require(0);
  196.         where.h=100;
  197.         where.v=100;
  198.         typeList[0]='TEXT';
  199.         reply.version=0;
  200.         SFGetFile(where,"\p",NULL,1,typeList,NULL,&reply);
  201.         if(!reply.good)PrintfExit("Couldn't open file.\n");
  202.         SetVol(NULL,reply.vRefNum);    /* look in that folder */
  203.         sprintf(dataFileName,"%#s",reply.fName);
  204.     #else
  205.         // Insert code here to get data file name into dataFileName.
  206.     #endif
  207.     dataFile=fopen(dataFileName,"r");
  208.     if(dataFile==NULL) {
  209.         PrintfExit("\007Sorry, can't find file \"%s\".\007\n",dataFileName);
  210.     }
  211.     printf("Reading \"%s\"\n",dataFileName);
  212.     strcpy(string,dataFileName);
  213.     if(strstr(string,".data"))strcpy(strstr(string,".data"),"");
  214.     else if(strstr(string,".yes"))strcpy(strstr(string,".yes"),"");
  215.     if(strlen(string)>0){
  216.         sprintf(fitFileName,"%s.fit",string);
  217.         printf("Creating output files:\n%s\n%s.<condition name>.plot\n",fitFileName,string);
  218.         fitFile=fopen(fitFileName,"w");
  219.         if(fitFile==NULL)PrintfExit("Sorry, I can't create file \"%s\".\007\n",fitFileName);
  220.         #if MACINTOSH
  221.             SetFileInfo(fitFileName,'TEXT','XCEL');        /* Excel file */
  222.         #endif
  223.     }
  224.     modelDF=PARAMS;    /* number of parameters to be adjusted in fitting */
  225.     printf("How many parameters shall be adjustable? 0 to %d (%d):",PARAMS,modelDF);
  226.     gets(string);
  227.     sscanf(string,"%d",&modelDF);
  228.     for(i=1;i<modelDF;i++){
  229.         printf("Initial value for %s? (%.2f):",paramName[i],paramPtr[i]);
  230.         gets(string);
  231.         sscanf(string,"%lf",¶mPtr[i]);
  232.     }
  233.     for(i=modelDF;i<PARAMS;i++){
  234.         printf("Fixed value for %s? (%.2f):",paramName[i],paramPtr[i]);
  235.         gets(string);
  236.         sscanf(string,"%lf",¶mPtr[i]);
  237.     }
  238.     initParams=params;
  239.     
  240.     printf("Reading %s\n",dataFileName);
  241.     if(fitFile)fprintf(fitFile,"Condition\t"
  242.         "%s\t%s\t%s\t%s\tfree params"
  243.         "\tsignif.\tChi sq.\td.f.\ttrials\tcontrasts\n",
  244.         paramName[0],paramName[1],paramName[2],paramName[3]);
  245.     while(!feof(dataFile)){
  246.         MyFgets(string,sizeof(string),dataFile);
  247.         printf("\n");
  248.         while (a=fgetc(dataFile),ungetc(a,dataFile),a=='#' || a==EOF){
  249.             printf("%s\n",string);                                /* echo comment */
  250.             if(fitFile)fprintf(fitFile,"%s\n",string);
  251.             MyFgets(string,sizeof(string),dataFile);
  252.             if(feof(dataFile))break;
  253.         }
  254.         /* get condition name. Strip leading # and trailing nonprinting junk */
  255.         strcpy(conditionName,&string[1]);
  256.         for(i=0;i<strlen(conditionName);i++){
  257.             if(!isprint(conditionName[i]))conditionName[i]=0;
  258.         }
  259.         printf("\n");
  260.         data.contrasts=i=0;
  261.         while(a=fgetc(dataFile),ungetc(a,dataFile),a!='#' && a!=EOF){    /* read data */
  262.             if(i==MAX_CONTRASTS){
  263.                 SortAndMergeContrasts(&data);
  264.                 if(i<MAX_CONTRASTS)break;
  265.                 PrintfExit("Quick3: Sorry, too many different contrasts. >MAX_CONTRASTS==%d\007\n",
  266.                     MAX_CONTRASTS);
  267.             }
  268.             MyFgets(string,sizeof(string),dataFile);
  269.             sscanf(string,"%lf\t%ld\t%ld",
  270.                 &data.c[i].contrast,&data.c[i].trials,&data.c[i].correct);
  271.             i++;
  272.             data.contrasts=i;
  273.         }
  274.         if(data.contrasts==0)break;
  275.         SortAndMergeContrasts(&data);
  276.         trials=0;
  277.         for(i=0;i<data.contrasts;i++)trials+=data.c[i].trials;    /* total trials */
  278.         params=initParams;    /* initial guess & fixed values */
  279.         if(modelDF>0)
  280.             params.logAlpha=log(data.c[data.contrasts/2].contrast)/log(10.0);    /* median contrast */
  281.         printf("%s\n",conditionName);
  282.         printf("\t%s   %s   %s   %s  contrasts trials signif. Chi sq.   d.f.\n",
  283.             paramName[0],paramName[1],paramName[2],paramName[3]);
  284.         printf("Guess:");
  285.         printf("\t%7.2f%8.1f%8.2f%8.2f\n",
  286.             paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3]);
  287.         significance=PsychometricFit(¶ms,ModelFunction,&data,&modelLL,modelDF,
  288.             &chiSquare,&chiSquareDF);
  289.         printf("Fit:");
  290.         printf("\t%7.2f%8.1f%8.2f%8.2f",
  291.             paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3]);
  292.         printf("%8d%8ld",data.contrasts,trials);
  293.         printf("%8.2f%8.1f%8d\n",significance,chiSquare,chiSquareDF);
  294.         if(fitFile)fprintf(fitFile,"%s"
  295.                 "\t%7.4f\t%7.4f\t%7.4f\t%7.4f\t%7d"
  296.                 "\t%7.4f\t%7.4f\t%7d"
  297.                 "\t%7ld\t%7d\n",
  298.                 conditionName,
  299.                 paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3],modelDF,
  300.                 significance,chiSquare,chiSquareDF,
  301.                 trials,data.contrasts);
  302.     
  303.         /* Now create plot file */
  304.         if(fitFile){
  305.             strcpy(plotFileName,fitFileName);
  306.             plotFileName[strlen(plotFileName)-strlen(".fit")]=0;    /* strip off the ".fit" */
  307.             sprintf(plotFileName,"%s.%s.plot",plotFileName,conditionName);
  308.             plotFile=fopen(plotFileName,"w");
  309.             if(plotFile==NULL){
  310.                 PrintfExit("Sorry, I can't create file \"%s\".\n\007",plotFileName);
  311.             }
  312.             #if MACINTOSH
  313.                 SetFileInfo(plotFileName,'TEXT','CGRF');    /* Cricket Graph file */
  314.             #endif
  315.             monotonicData=data;
  316.             MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  317.             fprintf(plotFile,"*\n");
  318.             fprintf(plotFile,"Contrast\t%s\tLower bound\tUpper bound\t%s fit\tMonotone fit\tTrials\tCorrect\tParameters\n",
  319.                 conditionName,modelName);
  320.             for(i=0;i<data.contrasts;i++){
  321.                 cPtr=&data.c[i];
  322.                 fprintf(plotFile,"%.3f",cPtr->contrast);
  323.                 if(cPtr->trials>0){
  324.                     fprintf(plotFile,"\t%.3f",cPtr->correct/(double)cPtr->trials);
  325.                     fprintf(plotFile,"\t%.3f\t%.3f",
  326.                         BinomialLowerBound(0.95,cPtr->correct,cPtr->trials),
  327.                         BinomialUpperBound(0.95,cPtr->correct,cPtr->trials));
  328.                 }
  329.                 else
  330.                     fprintf(plotFile,"\t\t\t");
  331.                 fprintf(plotFile,"\t%.3f",(*ModelFunction)(cPtr->contrast,¶ms));
  332.                 if(cPtr->trials>0)
  333.                     fprintf(plotFile,"\t%.3f\t%5ld\t%5ld",
  334.                         monotonicData.c[i].correct/(double)monotonicData.c[i].trials,
  335.                         cPtr->trials,cPtr->correct);
  336.                 else fprintf(plotFile,"\t\t\t");
  337.                 if(i<PARAMS)fprintf(plotFile,"\t%s=%.3f",paramName[i],paramPtr[i]);
  338.                 if(i==0)fprintf(plotFile,"\talpha=%.4f",pow(10.0,paramPtr[0]));
  339.                 if(i==PARAMS)fprintf(plotFile,"\tsignificance=%.3f",significance);
  340.                 fprintf(plotFile,"\n");
  341.             }
  342.             if(plotFile)fclose(plotFile);
  343.         }
  344.     }
  345.     if(fitFile)fclose(fitFile);
  346.     if(dataFile)fclose(dataFile);
  347. }
  348.  
  349.